
/*
 * Este programa calcula 3 autovalores de una matriz cuadrada:
 *  -> El autovalor de menor valor absoluto y su autovector asociado.
 *  -> El autovalor de mayor valor absoluto y su autovector asociado.
 *  -> Un autovalor de valor absoluto intermedio y con su autovector asociado.
 * 
 * La matriz de la cual se obtienen los autovalores puede obtenerse de 3 maneras:
 *  -> Especificando el nombre del archivo que contiene sus valores.
 *  -> Ingresando los valores manualmente al inicio del programa.
 *  -> Pidiendo al programa que genere una matriz aleatoria.
 *
 * Este programa realiza los cálculos en forma paralela mediante la librería mpich2.
 *
 * Para el cálculo de los autovalores y autovectores, se utiliza el método
 * de las potencias.
 *
 */

/* Inclusión de bibliotecas propias */
#include "matrices.h"
#include "autovalores.h"
#include "inversa.h"
#include "entradasalida.h"

/* Inclusión de bicliotecas estandares */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpi.h>

#define MAX_ITER 20000000
#define MIN_ERR 0.0000001


/*
 * Crea un objeto VV_PROPIO con un vector de dimensión especificada
 * y con autovalor 0.
 */

struct VV_PROPIO *crear_VV_PROPIO(int dimension);
int analizarTriangularizacion(struct MATRIZ *matriz);
struct MATRIZ *invertirMatrizL(struct MATRIZ *matrizL);
struct MATRIZ *invertirMatrizU(struct MATRIZ *matrizU);


int sustitucionInversaParalela(struct MATRIZ *coefs, double *ti, double *solucion) {

    // Buffereamos la dimensión de la matriz
    int dim = coefs->dimension;

    int i,j; // Variables de iteración
    double temp; // Variable buffer

    // Si el coeficiente de la incógnita es 0, no se puede determinar la inversa
    if (coefs->elementos[dim-1][dim-1] == 0) {
        return 0;
    }

    // Obtenemos la última incógnita
    solucion[dim-1] = ti[dim-1] / coefs->elementos[dim-1][dim-1];

    // Resolvemos las siguientes incógnitas
    for (i=dim-2; i>=0; i--) {

        // Si el coeficiente de la incógnita es 0, no se puede determinar la inversa
        if(coefs->elementos[i][i] == 0) {
           return 0;
        }

        temp = 0; // Reseteamos el buffer

        // Calculamos el valor de la incógnita actual
        for (j=i+1; j<dim; j++) {
            temp += coefs->elementos[i][j] * solucion[j];
        }

        // Asignamos el valor encontrado
        solucion[i] = (ti[i] - temp ) / coefs->elementos[i][i];

    }

    // Pudimos calcular todas las incógnitas
    return 1;

}


int sustitucionDirectaParalela(struct MATRIZ *coefs, double *ti, double *solucion) {

    int i,j; // Variables de iteración
    double temp; // Buffer

    // Si el primer elemento de la matriz de coeficientes es 0,
    // no se puede determinar la inversa.
    if(coefs->elementos[0][0] == 0) {
        return 0;
    }

    // Obtenemos la primera incognita.
    solucion[0] = ti[0] / coefs->elementos[0][0];

    // Resolver las siguients incógnitas.
    for (i=1; i < coefs->dimension; i++) {

        // Si el coeficiente de la incógnita es 0, no se puede determinar la inversa
        if(coefs->elementos[i][i] == 0) {
           return 0;
        }

        temp = 0; // Reiniciamos el buffer

        // Calculamos el valor de la incógnita actual
        for (j=0; j<i; j++) {
            temp += coefs->elementos[i][j] * solucion[j];
        }

        // Asignamos la incónita calculada
        solucion[i] = (ti[i] - temp ) / coefs->elementos[i][i];

    }

    // Pudimos calcular todas las incógnitas
    return 1;

}



/* Sección de implementación de funciones */

int main(int argc, char **argv) {

    // Iniciamos la librería mpich2
    MPI_Init(&argc,&argv);

    struct MATRIZ *matriz; // Aqui guardamos la matriz sobre la que operaremos
    time_t instanteInicial; // Instante inicial del calculo
    time_t instanteFinal; // Instante final del calculo
    struct VV_PROPIO *mayorVV;
    struct VV_PROPIO *menorVV;
    struct VV_PROPIO *medioVV;
    int globalRank;

    // Verificar nuestro ranking dentro del comunicador global
    MPI_Comm_rank(MPI_COMM_WORLD, &globalRank);

    if(globalRank == 0) {
        
        // Imprimimos la bienvenida
        printf("============================================================================\n");
        printf("Cálculo del menor/medio/mayor autovalor y su autovector. VERSIÓN PARALELA\n");
        printf("============================================================================\n");

        // Aquí guardaremos la matriz con la cual haremos los cálculos y su tamaño.
        matriz = obtenerMatrizDelUsuario();

        // Guardamos el momento inicial
        instanteInicial = time(NULL);

    }

    /*
     * PRIMERA PARTE DEL PROCESO PARALELO:
     * Calculo del mayor autovalor
     */

     // Declaración de variables
    int tamanio; // Tamaño de las matrices (n x n)
    long contadorIteraciones;  // Contador de la cantidad de iteraciones que lleva el método
    double lambdaActual;      // Buffer del autovalor obtenido en la iteracion actual
    double lambdaAnterior;  // Buffer del autovalor obtenido en la iteracion anterior a la actual
    struct VECTOR *resultado; // Buffer del calculo de A*b
    int pivote;

    // Creamos el objecto vv_propio y obtenemos su vector interno
    struct VV_PROPIO *vvPropio;
    struct VECTOR *autoVector;
    int i,j;
    int *cantidadFilas;
    int condicionSalida; // Aquí guardaremos si se debe terminar de iterar o no.
    double errorActual; // VAriación del autovalor con respecto a la iteración anterior

    int cantidadWorkers;

    int cantidadFilasNodo;
    double **filasNodo;
    double *resultadoNodo;

    if (globalRank==0) {

        tamanio = matriz->dimension;
        contadorIteraciones = 0;
        lambdaActual = 0;
        lambdaAnterior = 0;
        resultado = crearVector(tamanio);
        pivote = -1;
        vvPropio = crear_VV_PROPIO(tamanio);
        autoVector = vvPropio->autovector;

        // Inicializamos la autovector inicial como un vector de 1
        for (i=0;i<tamanio;i++) {
            autoVector->elementos[i] = 1;
        }

        // Decidir cuantas filas le enviamos a cada nodo

        // Vemos cuantos nodos hay y cuantos workers tenemos
        int cantidadNodos;
        MPI_Comm_size(MPI_COMM_WORLD, &cantidadNodos);
        cantidadWorkers = cantidadNodos -1;

        int cantidadInicial = tamanio / cantidadWorkers;
        int cantidadExtra = tamanio % cantidadWorkers;
        cantidadFilas = malloc(sizeof(int) * cantidadNodos);

        for(i=0; i<cantidadWorkers; i++) {
            cantidadFilas[i+1] = cantidadInicial;
        }

        for(i=0; i<cantidadExtra; i++) {
            cantidadFilas[i+1] += 1;
        }

    }

    // Difundimos el tamaño del nodo
    MPI_Bcast(&tamanio, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(globalRank != 0) {
        autoVector = crearVector(tamanio);
    }

    // Le decimos a cada nodo cuantas filas va a recibir
    MPI_Scatter(cantidadFilas, 1, MPI_INT, &cantidadFilasNodo, 1, MPI_INT, 0, MPI_COMM_WORLD);

    // Enviamos las filas a cada nodo
    if (globalRank == 0) {

        int contadorFilas = 0;

        for(i=0; i<cantidadWorkers; i++) {
            
            for (j=0; j<cantidadFilas[i+1]; j++) {
                MPI_Send(matriz->elementos[contadorFilas], tamanio, MPI_DOUBLE, i+1, 1000, MPI_COMM_WORLD);
                contadorFilas++;
            }
            
        }
        
    } else {

        filasNodo = malloc(sizeof(double *) * cantidadFilasNodo);
        resultadoNodo = malloc(sizeof(double) * cantidadFilasNodo);

        for(i=0; i<cantidadFilasNodo; i++) {
            filasNodo[i] = malloc(sizeof(double) * tamanio);
            MPI_Status status;
            MPI_Recv(filasNodo[i], tamanio, MPI_DOUBLE, 0, 1000, MPI_COMM_WORLD, &status);
        }

    }

    do {

        if (globalRank == 0) {

            // Buffereamos el lambda de la iteracion anterior
            lambdaAnterior = lambdaActual;

        }

        // Difundimos el vector de la iteracion actual
        MPI_Bcast(autoVector->elementos, tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        if(globalRank == 0) {

            int resultadoOffset = 0;

            // ESPERAR EL RESULTADO
            for(i=0; i< cantidadWorkers; i++) {

                if(cantidadFilas[i+1] != 0) {
                    MPI_Status status;
                    MPI_Recv(resultado->elementos + resultadoOffset, cantidadFilas[i+1], MPI_DOUBLE, i+1, 2000+i+1, MPI_COMM_WORLD, &status);
                    resultadoOffset += cantidadFilas[i+1];
                }

            }

        } else {

            int k;
            
            // MULTIPLICAR
            for (k=0; k<cantidadFilasNodo; k++) {

                resultadoNodo[k] = 0;

                for (j=0; j<tamanio; j++) {
                    resultadoNodo[k] += filasNodo[k][j] * autoVector->elementos[j];
                }

            }

            if (cantidadFilasNodo != 0) {

                // Enviamos los resultados
                MPI_Send(resultadoNodo, cantidadFilasNodo, MPI_DOUBLE, 0, 2000+globalRank, MPI_COMM_WORLD);

            }

        }

        if(globalRank == 0) {

            // Seleccionamos un elemento como pivot (si no se hizo previamente)
            if (pivote<0) {
                pivote = buscarMayorElemento(resultado);
            }
            
            // Seleccionamos el nuevo autovalor a partir del pivote seleccionado
            lambdaActual = resultado->elementos[pivote];

            // Calculamos el nuevo autovector
            for (i=0;i<tamanio;i++) {
                autoVector->elementos[i] = resultado->elementos[i] / lambdaActual;
            }

            // Calculamos el error actual
            errorActual = fabs((lambdaActual - lambdaAnterior) / lambdaActual);

            // Incrementamos la cantidad de iteraciones en 1
            contadorIteraciones++;

            /*
             * Iteramos hasta que se cumpla una de la siguientes condiciones:
             * 1) La diferencia porcentual entre el autovalor anterior y el actual
             *    sea menor que el error indicado por parámetro.
             * 2) Alcancemos el límite de iteraciones especificado por parámetro.
             */
            condicionSalida = errorActual <= MIN_ERR
                    || contadorIteraciones == MAX_ITER;

        }

        MPI_Bcast(&condicionSalida, 1, MPI_INT, 0, MPI_COMM_WORLD);

    } while (!condicionSalida);

    if(globalRank == 0) {

        // Liberamos la memoria reservada
        eliminarVector(resultado);

        // Seteamos el autovalor calculado en el objeto VV_PROPIO
        vvPropio->autovalor = lambdaActual;

        mayorVV = vvPropio;

    }

/******************************** FIN DE LA PRIMERA PARTE ************************************/

    MPI_Barrier(MPI_COMM_WORLD);

/******************************** COMIENZO DE LA SEGUNDA PARTE *******************************/

    struct MATRIZ *matrizL;
    struct MATRIZ *matrizU;
    struct MATRIZ *inversaL;
    struct MATRIZ *inversaU;
    struct MATRIZ *matrizInvertida;

    if (globalRank == 0) {

        // Si no se puede triangularizar, devolvemos NULL
        if(!analizarTriangularizacion(matriz)) {
            printf("\nERROR: La matriz no se puede triangularizar");
            MPI_Finalize();
            return MPI_SUCCESS;
        }

        // Creamos las matrices L y U.
        matrizL = crearMatriz(matriz->dimension);
        matrizU = crearMatriz(matriz->dimension);

        // Inicializamos ambas matrices
        int i,j;
        for (i=0; i < tamanio; i++) {

            // Inicializamos U
            for (j=0; j < tamanio ; j++) {
                matrizU->elementos[i][j] = (matriz->elementos)[i][j];
            }

            // Inicializamos L
            matrizL->elementos[i][i] = 1;
            for (j=i+1; j<tamanio; j++) {
                matrizL->elementos[i][j] = 0;
            }

        }

    }

    double multiplicador; // Buffer del multiplicador por fila
    int k,w,h; // Iterador por elementos
    int *filasNodoO;
    int filasTriangNodo;
    double *filaPivote;
    double **filasT;

    if (globalRank == 0) {
        filasNodoO = malloc(sizeof(int) * cantidadWorkers+1);
    }

    // Recorremos la matriz por columnas
    for (i=0; i<tamanio-1; i++) {

        if(globalRank == 0) {

            // Decidimos cuantas filas le enviaremos a cada nodo
            int cantidadInicial = (tamanio-i-1) / cantidadWorkers;
            int cantidadExtra = (tamanio-i-1) % cantidadWorkers;

            for(w=0; w<cantidadWorkers; w++) {
                filasNodoO[w+1] = cantidadInicial;
            }

            for(w=0; w<cantidadExtra; w++) {
                filasNodoO[w+1] += 1;
            }

            filaPivote = matrizU->elementos[i];

        } else {
            filaPivote = malloc(sizeof(double) * tamanio);
        }

        // Avisamos cuantas filas enviaremos a cada nodo
        MPI_Scatter(filasNodoO, 1, MPI_INT, &filasTriangNodo, 1, MPI_INT, 0, MPI_COMM_WORLD);

        // Distribuimos la fila pivote
        MPI_Bcast(filaPivote, tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        if(globalRank == 0) {

            int offsetFilas = i+1;

            // Distribuimos las filas correspondientes a cada nodo
            for(w=0; w<cantidadWorkers; w++) {
                for(j=0; j<filasNodoO[w+1] ;j++) {
                    MPI_Send(matrizU->elementos[offsetFilas], tamanio, MPI_DOUBLE, w+1, 3001+w, MPI_COMM_WORLD);
                    offsetFilas += 1;
                }
            }

        } else {

            // Recibimos las filas correspondientes
            filasT = malloc(sizeof(double *) * filasTriangNodo);
            for(j=0; j<filasTriangNodo; j++) {
                filasT[j] = malloc(sizeof(double) * tamanio);
                MPI_Status status;
                MPI_Recv(filasT[j], tamanio, MPI_DOUBLE, 0, 3000+globalRank, MPI_COMM_WORLD, &status);
            }
            
        }

        if (globalRank == 0) {

            int offsetMatriz = i+1;

            // Recibimos los resultados de los nodos
            for(j=0; j<cantidadWorkers; j++) {
                for(h=0; h<filasNodoO[j+1]; h++) {
                    MPI_Status status;
                    MPI_Recv(matrizU->elementos[offsetMatriz], tamanio, MPI_DOUBLE, j+1, 4000+h, MPI_COMM_WORLD, &status);
                    offsetMatriz += 1;
                }
            }

        } else {

            for (j=0; j<filasTriangNodo; j++) {

                // Calculo el multiplicador del elemento i,j.
                multiplicador = filasT[j][i] / filaPivote[i];

                // Modificamos la fila j
                for (k=i; k<tamanio; k++){
                    filasT[j][k] -= multiplicador * filaPivote[k];
                }

                // Guardamos el multiplicador
                filasT[j][i] = multiplicador;

            }

            // Enviamos los resultados al master
            for(j=0; j<filasTriangNodo; j++) {
                MPI_Send(filasT[j], tamanio, MPI_DOUBLE, 0, 4000+j, MPI_COMM_WORLD);
            }

        }

    }

    if(globalRank == 0) {

        // Extraemos la matrices L y U
        for (i=1; i<tamanio; i++) {
            for (j=0; j<i; j++) {
                matrizL->elementos[i][j] = matrizU->elementos[i][j];
                matrizU->elementos[i][j] = 0;
            }
        }

    }

    // ACA COMIENZA EL CALCULO DE LA MATRIZ INVERSA

    if (globalRank == 0) {

        // Aquí guardaremos la inversa de la matriz
        inversaL = crearMatriz(tamanio);

    } else {
        matrizL = crearMatriz(tamanio);
    }

    // Difundir la matriz L
    for(i=0; i<tamanio; i++) {
        MPI_Bcast(matrizL->elementos[i], tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    }

    int *indicesI;

    // Nuevamente decidimos cuantas filas enviar a cada proceso
    if (globalRank == 0) {

        indicesI = malloc(sizeof(int) * (cantidadWorkers + 1));
        int contadorFil = 0;

        // Decidimos cuantas filas le enviaremos a cada nodo
        int cantidadInicial = (tamanio) / cantidadWorkers;
        int cantidadExtra = (tamanio) % cantidadWorkers;

        for(w=0; w<cantidadWorkers; w++) {
            filasNodoO[w+1] = cantidadInicial;
        }

        for(w=0; w<cantidadExtra; w++) {
            filasNodoO[w+1] += 1;
        }

        for(w=0; w<cantidadWorkers; w++) {
            indicesI[w+1] = contadorFil;
            contadorFil += filasNodoO[w+1];
        }

    }

    int indiceInicial;
    MPI_Scatter(indicesI, 1, MPI_INT, &indiceInicial, 1, MPI_INT, 0, MPI_COMM_WORLD);

    int cuantasFilasU;
    MPI_Scatter(filasNodoO, 1, MPI_INT, &cuantasFilasU, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(globalRank == 0) {

        int offsetMatriz = 0;

        double *queBuffer = malloc(sizeof(double) * tamanio);

        // Recibimos los resultados de los nodos
        for(j=0; j<cantidadWorkers; j++) {
            for(h=0; h<filasNodoO[j+1]; h++) {
                MPI_Status status;
                MPI_Recv(queBuffer, tamanio, MPI_DOUBLE, j+1, 5000+h, MPI_COMM_WORLD, &status);
                for(w=0; w<tamanio; w++) {
                    inversaL->elementos[w][offsetMatriz] = queBuffer[w];
                }
                offsetMatriz += 1;
            }
        }

    } else {

        // El worker hace su trabajo
        int columnaInicial = indiceInicial;
        double **misColumnas = malloc(sizeof(double *) * cuantasFilasU);
        double **misRespuestas = malloc(sizeof(double *) * cuantasFilasU);
        
        for (w=0; w<cuantasFilasU; w++) {
            misColumnas[w] = malloc(sizeof(double) * tamanio);
            misRespuestas[w] = malloc(sizeof(double) * tamanio);
            for (h=0; h<tamanio; h++) {
                if (h==columnaInicial) {
                    misColumnas[w][h] = 1;
                } else {
                    misColumnas[w][h] = 0;
                }
            }
            columnaInicial++;
        }

        for (w=0; w<cuantasFilasU; w++) {
            sustitucionDirectaParalela(matrizL,misColumnas[w],misRespuestas[w]);
        }

        for (w=0; w<cuantasFilasU; w++) {
            MPI_Send(misRespuestas[w], tamanio, MPI_DOUBLE, 0, 5000+w, MPI_COMM_WORLD);
        }

    }


    // ACA COMIENZA EL CALCULO DE LA MATRIZ INVERSA

    if (globalRank == 0) {

        // Aquí guardaremos la inversa de la matriz
        inversaU = crearMatriz(tamanio);

    } else {
        matrizU = crearMatriz(tamanio);
    }

    // Difundir la matriz L
    for(i=0; i<tamanio; i++) {
        MPI_Bcast(matrizU->elementos[i], tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    }

    // Nuevamente decidimos cuantas filas enviar a cada proceso
    if (globalRank == 0) {

        int contadorFil = 0;

        // Decidimos cuantas filas le enviaremos a cada nodo
        int cantidadInicial = (tamanio) / cantidadWorkers;
        int cantidadExtra = (tamanio) % cantidadWorkers;

        for(w=0; w<cantidadWorkers; w++) {
            filasNodoO[w+1] = cantidadInicial;
        }

        for(w=0; w<cantidadExtra; w++) {
            filasNodoO[w+1] += 1;
        }

        for(w=0; w<cantidadWorkers; w++) {
            indicesI[w+1] = contadorFil;
            contadorFil += filasNodoO[w+1];
        }

    }

    MPI_Scatter(indicesI, 1, MPI_INT, &indiceInicial, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Scatter(filasNodoO, 1, MPI_INT, &cuantasFilasU, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(globalRank == 0) {

        int offsetMatriz = 0;

        double *queBuffer = malloc(sizeof(double) * tamanio);

        // Recibimos los resultados de los nodos
        for(j=0; j<cantidadWorkers; j++) {
            for(h=0; h<filasNodoO[j+1]; h++) {
                MPI_Status status;
                MPI_Recv(queBuffer, tamanio, MPI_DOUBLE, j+1, 6000+h, MPI_COMM_WORLD, &status);
                for(w=0; w<tamanio; w++) {
                    inversaU->elementos[w][offsetMatriz] = queBuffer[w];
                }
                offsetMatriz += 1;
            }
        }

    } else {

        // El worker hace su trabajo
        int columnaInicial = indiceInicial;
        double **misColumnas = malloc(sizeof(double *) * cuantasFilasU);
        double **misRespuestas = malloc(sizeof(double *) * cuantasFilasU);

        for (w=0; w<cuantasFilasU; w++) {
            misColumnas[w] = malloc(sizeof(double) * tamanio);
            misRespuestas[w] = malloc(sizeof(double) * tamanio);
            for (h=0; h<tamanio; h++) {
                if (h==columnaInicial) {
                    misColumnas[w][h] = 1;
                } else {
                    misColumnas[w][h] = 0;
                }
            }
            columnaInicial++;
        }

        for (w=0; w<cuantasFilasU; w++) {
            sustitucionInversaParalela(matrizU,misColumnas[w],misRespuestas[w]);
        }

        for (w=0; w<cuantasFilasU; w++) {
            MPI_Send(misRespuestas[w], tamanio, MPI_DOUBLE, 0, 6000+w, MPI_COMM_WORLD);
        }

    }    


    if(globalRank==0) {

        // Si algunas de las matrices es nula, devolver NULL
        if(inversaL == NULL || inversaU == NULL) {
            printf("\nNo se pudo calcular la inversa de la matriz");
            MPI_Finalize();
            return MPI_SUCCESS;
        }

        // Multiplicamos ambas inversas y obtenemos la inversa de la matriz
        matrizInvertida = multiplicarMatricesLU(inversaU,inversaL);

        // Eliminamos las inversas de L y U
        eliminarMatriz(inversaL);
        eliminarMatriz(inversaU);

    }

/******************************** FIN DE LA SEGUNDA PARTE ************************************/

    MPI_Barrier(MPI_COMM_WORLD);

/******************************** COMIENZO DE LA TERCERA PARTE *******************************/

    if (globalRank == 0) {

        // Si la matriz devuelta fue NULL, salir
        if (matrizInvertida==NULL) {
            printf("\nNo se pudo calcular la inversa de la matriz principal.\n\n");
            return EXIT_FAILURE;
        }

    }

    if (globalRank==0) {

        contadorIteraciones = 0;
        lambdaActual = 0;
        lambdaAnterior = 0;
        resultado = crearVector(tamanio);
        pivote = -1;
        vvPropio = crear_VV_PROPIO(tamanio);
        autoVector = vvPropio->autovector;

        // Inicializamos la autovector inicial como un vector de 1
        for (i=0;i<tamanio;i++) {
            autoVector->elementos[i] = 1;
        }

        // Decidir cuantas filas le enviamos a cada nodo

        // Vemos cuantos nodos hay y cuantos workers tenemos
        int cantidadNodos;
        MPI_Comm_size(MPI_COMM_WORLD, &cantidadNodos);
        cantidadWorkers = cantidadNodos -1;

        int cantidadInicial = tamanio / cantidadWorkers;
        int cantidadExtra = tamanio % cantidadWorkers;
        cantidadFilas = malloc(sizeof(int) * cantidadNodos);

        for(i=0; i<cantidadWorkers; i++) {
            cantidadFilas[i+1] = cantidadInicial;
        }

        for(i=0; i<cantidadExtra; i++) {
            cantidadFilas[i+1] += 1;
        }

    }

    // Difundimos el tamaño del nodo
    MPI_Bcast(&tamanio, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(globalRank != 0) {
        autoVector = crearVector(tamanio);
    }

    // Le decimos a cada nodo cuantas filas va a recibir
    MPI_Scatter(cantidadFilas, 1, MPI_INT, &cantidadFilasNodo, 1, MPI_INT, 0, MPI_COMM_WORLD);

    // Enviamos las filas a cada nodo
    if (globalRank == 0) {

        int contadorFilas = 0;

        for(i=0; i<cantidadWorkers; i++) {

            for (j=0; j<cantidadFilas[i+1]; j++) {
                MPI_Send(matrizInvertida->elementos[contadorFilas], tamanio, MPI_DOUBLE, i+1, 1000, MPI_COMM_WORLD);
                contadorFilas++;
            }

        }

    } else {

        filasNodo = malloc(sizeof(double *) * cantidadFilasNodo);
        resultadoNodo = malloc(sizeof(double) * cantidadFilasNodo);

        for(i=0; i<cantidadFilasNodo; i++) {
            filasNodo[i] = malloc(sizeof(double) * tamanio);
            MPI_Status status;
            MPI_Recv(filasNodo[i], tamanio, MPI_DOUBLE, 0, 1000, MPI_COMM_WORLD, &status);
        }

    }

    do {

        if (globalRank == 0) {

            // Buffereamos el lambda de la iteracion anterior
            lambdaAnterior = lambdaActual;

        }

        // Difundimos el vector de la iteracion actual
        MPI_Bcast(autoVector->elementos, tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        if(globalRank == 0) {

            int resultadoOffset = 0;

            // ESPERAR EL RESULTADO
            for(i=0; i< cantidadWorkers; i++) {

                if(cantidadFilas[i+1] != 0) {
                    MPI_Status status;
                    MPI_Recv(resultado->elementos + resultadoOffset, cantidadFilas[i+1], MPI_DOUBLE, i+1, 2000+i+1, MPI_COMM_WORLD, &status);
                    resultadoOffset += cantidadFilas[i+1];
                }

            }

        } else {

            int k;

            // MULTIPLICAR
            for (k=0; k<cantidadFilasNodo; k++) {

                resultadoNodo[k] = 0;

                for (j=0; j<tamanio; j++) {
                    resultadoNodo[k] += filasNodo[k][j] * autoVector->elementos[j];
                }

            }

            if (cantidadFilasNodo != 0) {

                // Enviamos los resultados
                MPI_Send(resultadoNodo, cantidadFilasNodo, MPI_DOUBLE, 0, 2000+globalRank, MPI_COMM_WORLD);

            }

        }

        if(globalRank == 0) {

            // Seleccionamos un elemento como pivot (si no se hizo previamente)
            if (pivote<0) {
                pivote = buscarMayorElemento(resultado);
            }

            // Seleccionamos el nuevo autovalor a partir del pivote seleccionado
            lambdaActual = resultado->elementos[pivote];

            // Calculamos el nuevo autovector
            for (i=0;i<tamanio;i++) {
                autoVector->elementos[i] = resultado->elementos[i] / lambdaActual;
            }

            // Calculamos el error actual
            errorActual = fabs((lambdaActual - lambdaAnterior) / lambdaActual);

            // Incrementamos la cantidad de iteraciones en 1
            contadorIteraciones++;

            /*
             * Iteramos hasta que se cumpla una de la siguientes condiciones:
             * 1) La diferencia porcentual entre el autovalor anterior y el actual
             *    sea menor que el error indicado por parámetro.
             * 2) Alcancemos el límite de iteraciones especificado por parámetro.
             */
            condicionSalida = errorActual <= MIN_ERR
                    || contadorIteraciones == MAX_ITER;

        }

        MPI_Bcast(&condicionSalida, 1, MPI_INT, 0, MPI_COMM_WORLD);

    } while (!condicionSalida);

    if(globalRank == 0) {

        // Liberamos la memoria reservada
        eliminarVector(resultado);

        // Seteamos el autovalor calculado en el objeto VV_PROPIO
        vvPropio->autovalor = lambdaActual;

        menorVV = vvPropio;
        menorVV->autovalor = 1 / menorVV->autovalor;

        // Ya no necesitamos la matriz invertida
        eliminarMatriz(matrizInvertida);

    }


    MPI_Barrier(MPI_COMM_WORLD);

    struct MATRIZ *matrizModificada;
    double s;
    
    if (globalRank == 0) {

        // Obtenemos un valor medio entre los autovalores extremos
        s = aproximarAutovalorMedio(mayorVV->autovalor,menorVV->autovalor);

        // Creamos y calculamos la matriz modificada
        printf("Calculamos la matriz para el autovalor medio\n");
        restarIxs(matriz,s);

    }


/********************************************************************************************************/

    if (globalRank == 0) {

        // Si no se puede triangularizar, devolvemos NULL
        if(!analizarTriangularizacion(matriz)) {
            printf("\nERROR: La matriz no se puede triangularizar");
            MPI_Finalize();
            return MPI_SUCCESS;
        }

        // Creamos las matrices L y U.
        matrizL = crearMatriz(matriz->dimension);
        matrizU = crearMatriz(matriz->dimension);

        // Inicializamos ambas matrices
        int i,j;
        for (i=0; i < tamanio; i++) {

            // Inicializamos U
            for (j=0; j < tamanio ; j++) {
                matrizU->elementos[i][j] = (matriz->elementos)[i][j];
            }

            // Inicializamos L
            matrizL->elementos[i][i] = 1;
            for (j=i+1; j<tamanio; j++) {
                matrizL->elementos[i][j] = 0;
            }

        }

    }

    if (globalRank == 0) {
        filasNodoO = malloc(sizeof(int) * cantidadWorkers+1);
    }

    // Recorremos la matriz por columnas
    for (i=0; i<tamanio-1; i++) {

        if(globalRank == 0) {

            // Decidimos cuantas filas le enviaremos a cada nodo
            int cantidadInicial = (tamanio-i-1) / cantidadWorkers;
            int cantidadExtra = (tamanio-i-1) % cantidadWorkers;

            for(w=0; w<cantidadWorkers; w++) {
                filasNodoO[w+1] = cantidadInicial;
            }

            for(w=0; w<cantidadExtra; w++) {
                filasNodoO[w+1] += 1;
            }

            filaPivote = matrizU->elementos[i];

        } else {
            filaPivote = malloc(sizeof(double) * tamanio);
        }

        // Avisamos cuantas filas enviaremos a cada nodo
        MPI_Scatter(filasNodoO, 1, MPI_INT, &filasTriangNodo, 1, MPI_INT, 0, MPI_COMM_WORLD);

        // Distribuimos la fila pivote
        MPI_Bcast(filaPivote, tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        if(globalRank == 0) {

            int offsetFilas = i+1;

            // Distribuimos las filas correspondientes a cada nodo
            for(w=0; w<cantidadWorkers; w++) {
                for(j=0; j<filasNodoO[w+1] ;j++) {
                    MPI_Send(matrizU->elementos[offsetFilas], tamanio, MPI_DOUBLE, w+1, 3001+w, MPI_COMM_WORLD);
                    offsetFilas += 1;
                }
            }

        } else {

            // Recibimos las filas correspondientes
            filasT = malloc(sizeof(double *) * filasTriangNodo);
            for(j=0; j<filasTriangNodo; j++) {
                filasT[j] = malloc(sizeof(double) * tamanio);
                MPI_Status status;
                MPI_Recv(filasT[j], tamanio, MPI_DOUBLE, 0, 3000+globalRank, MPI_COMM_WORLD, &status);
            }

        }

        if (globalRank == 0) {

            int offsetMatriz = i+1;

            // Recibimos los resultados de los nodos
            for(j=0; j<cantidadWorkers; j++) {
                for(h=0; h<filasNodoO[j+1]; h++) {
                    MPI_Status status;
                    MPI_Recv(matrizU->elementos[offsetMatriz], tamanio, MPI_DOUBLE, j+1, 4000+h, MPI_COMM_WORLD, &status);
                    offsetMatriz += 1;
                }
            }

        } else {

            for (j=0; j<filasTriangNodo; j++) {

                // Calculo el multiplicador del elemento i,j.
                multiplicador = filasT[j][i] / filaPivote[i];

                // Modificamos la fila j
                for (k=i; k<tamanio; k++){
                    filasT[j][k] -= multiplicador * filaPivote[k];
                }

                // Guardamos el multiplicador
                filasT[j][i] = multiplicador;

            }

            // Enviamos los resultados al master
            for(j=0; j<filasTriangNodo; j++) {
                MPI_Send(filasT[j], tamanio, MPI_DOUBLE, 0, 4000+j, MPI_COMM_WORLD);
            }

        }

    }

    if(globalRank == 0) {

        // Extraemos la matrices L y U
        for (i=1; i<tamanio; i++) {
            for (j=0; j<i; j++) {
                matrizL->elementos[i][j] = matrizU->elementos[i][j];
                matrizU->elementos[i][j] = 0;
            }
        }

    }

    // ACA COMIENZA EL CALCULO DE LA MATRIZ INVERSA

    if (globalRank == 0) {

        // Aquí guardaremos la inversa de la matriz
        inversaL = crearMatriz(tamanio);

    } else {
        matrizL = crearMatriz(tamanio);
    }

    // Difundir la matriz L
    for(i=0; i<tamanio; i++) {
        MPI_Bcast(matrizL->elementos[i], tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    }

    // Nuevamente decidimos cuantas filas enviar a cada proceso
    if (globalRank == 0) {

        indicesI = malloc(sizeof(int) * (cantidadWorkers + 1));
        int contadorFil = 0;

        // Decidimos cuantas filas le enviaremos a cada nodo
        int cantidadInicial = (tamanio) / cantidadWorkers;
        int cantidadExtra = (tamanio) % cantidadWorkers;

        for(w=0; w<cantidadWorkers; w++) {
            filasNodoO[w+1] = cantidadInicial;
        }

        for(w=0; w<cantidadExtra; w++) {
            filasNodoO[w+1] += 1;
        }

        for(w=0; w<cantidadWorkers; w++) {
            indicesI[w+1] = contadorFil;
            contadorFil += filasNodoO[w+1];
        }

    }

    MPI_Scatter(indicesI, 1, MPI_INT, &indiceInicial, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Scatter(filasNodoO, 1, MPI_INT, &cuantasFilasU, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(globalRank == 0) {

        int offsetMatriz = 0;

        double *queBuffer = malloc(sizeof(double) * tamanio);

        // Recibimos los resultados de los nodos
        for(j=0; j<cantidadWorkers; j++) {
            for(h=0; h<filasNodoO[j+1]; h++) {
                MPI_Status status;
                MPI_Recv(queBuffer, tamanio, MPI_DOUBLE, j+1, 5000+h, MPI_COMM_WORLD, &status);
                for(w=0; w<tamanio; w++) {
                    inversaL->elementos[w][offsetMatriz] = queBuffer[w];
                }
                offsetMatriz += 1;
            }
        }

    } else {

        // El worker hace su trabajo
        int columnaInicial = indiceInicial;
        double **misColumnas = malloc(sizeof(double *) * cuantasFilasU);
        double **misRespuestas = malloc(sizeof(double *) * cuantasFilasU);

        for (w=0; w<cuantasFilasU; w++) {
            misColumnas[w] = malloc(sizeof(double) * tamanio);
            misRespuestas[w] = malloc(sizeof(double) * tamanio);
            for (h=0; h<tamanio; h++) {
                if (h==columnaInicial) {
                    misColumnas[w][h] = 1;
                } else {
                    misColumnas[w][h] = 0;
                }
            }
            columnaInicial++;
        }

        for (w=0; w<cuantasFilasU; w++) {
            sustitucionDirectaParalela(matrizL,misColumnas[w],misRespuestas[w]);
        }

        for (w=0; w<cuantasFilasU; w++) {
            MPI_Send(misRespuestas[w], tamanio, MPI_DOUBLE, 0, 5000+w, MPI_COMM_WORLD);
        }

    }


    // ACA COMIENZA EL CALCULO DE LA MATRIZ INVERSA

    if (globalRank == 0) {

        // Aquí guardaremos la inversa de la matriz
        inversaU = crearMatriz(tamanio);

    } else {
        matrizU = crearMatriz(tamanio);
    }

    // Difundir la matriz L
    for(i=0; i<tamanio; i++) {
        MPI_Bcast(matrizU->elementos[i], tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    }

    // Nuevamente decidimos cuantas filas enviar a cada proceso
    if (globalRank == 0) {

        int contadorFil = 0;

        // Decidimos cuantas filas le enviaremos a cada nodo
        int cantidadInicial = (tamanio) / cantidadWorkers;
        int cantidadExtra = (tamanio) % cantidadWorkers;

        for(w=0; w<cantidadWorkers; w++) {
            filasNodoO[w+1] = cantidadInicial;
        }

        for(w=0; w<cantidadExtra; w++) {
            filasNodoO[w+1] += 1;
        }

        for(w=0; w<cantidadWorkers; w++) {
            indicesI[w+1] = contadorFil;
            contadorFil += filasNodoO[w+1];
        }

    }

    MPI_Scatter(indicesI, 1, MPI_INT, &indiceInicial, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Scatter(filasNodoO, 1, MPI_INT, &cuantasFilasU, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(globalRank == 0) {

        int offsetMatriz = 0;

        double *queBuffer = malloc(sizeof(double) * tamanio);

        // Recibimos los resultados de los nodos
        for(j=0; j<cantidadWorkers; j++) {
            for(h=0; h<filasNodoO[j+1]; h++) {
                MPI_Status status;
                MPI_Recv(queBuffer, tamanio, MPI_DOUBLE, j+1, 6000+h, MPI_COMM_WORLD, &status);
                for(w=0; w<tamanio; w++) {
                    inversaU->elementos[w][offsetMatriz] = queBuffer[w];
                }
                offsetMatriz += 1;
            }
        }

    } else {

        // El worker hace su trabajo
        int columnaInicial = indiceInicial;
        double **misColumnas = malloc(sizeof(double *) * cuantasFilasU);
        double **misRespuestas = malloc(sizeof(double *) * cuantasFilasU);

        for (w=0; w<cuantasFilasU; w++) {
            misColumnas[w] = malloc(sizeof(double) * tamanio);
            misRespuestas[w] = malloc(sizeof(double) * tamanio);
            for (h=0; h<tamanio; h++) {
                if (h==columnaInicial) {
                    misColumnas[w][h] = 1;
                } else {
                    misColumnas[w][h] = 0;
                }
            }
            columnaInicial++;
        }

        for (w=0; w<cuantasFilasU; w++) {
            sustitucionInversaParalela(matrizU,misColumnas[w],misRespuestas[w]);
        }

        for (w=0; w<cuantasFilasU; w++) {
            MPI_Send(misRespuestas[w], tamanio, MPI_DOUBLE, 0, 6000+w, MPI_COMM_WORLD);
        }

    }


    if(globalRank==0) {

        // Si algunas de las matrices es nula, devolver NULL
        if(inversaL == NULL || inversaU == NULL) {
            printf("\nNo se pudo calcular la inversa de la matriz");
            MPI_Finalize();
            return MPI_SUCCESS;
        }

        // Multiplicamos ambas inversas y obtenemos la inversa de la matriz
        matrizModificada = multiplicarMatricesLU(inversaU,inversaL);

        // Eliminamos las inversas de L y U
        eliminarMatriz(inversaL);
        eliminarMatriz(inversaU);


        // Si la matriz devuelta fue NULL, salir
        if (matrizModificada==NULL) {
            printf("\nNo se pudo calcular la inversa de la matriz modificada.\n\n");
            return EXIT_FAILURE;
        }

        // Ya no necesitamos la matriz modificada
        eliminarMatriz(matriz);

    }

    if (globalRank==0) {

        contadorIteraciones = 0;
        lambdaActual = 0;
        lambdaAnterior = 0;
        resultado = crearVector(tamanio);
        pivote = -1;
        vvPropio = crear_VV_PROPIO(tamanio);
        autoVector = vvPropio->autovector;

        // Inicializamos la autovector inicial como un vector de 1
        for (i=0;i<tamanio;i++) {
            autoVector->elementos[i] = 1;
        }

        // Decidir cuantas filas le enviamos a cada nodo

        // Vemos cuantos nodos hay y cuantos workers tenemos
        int cantidadNodos;
        MPI_Comm_size(MPI_COMM_WORLD, &cantidadNodos);
        cantidadWorkers = cantidadNodos -1;

        int cantidadInicial = tamanio / cantidadWorkers;
        int cantidadExtra = tamanio % cantidadWorkers;
        cantidadFilas = malloc(sizeof(int) * cantidadNodos);

        for(i=0; i<cantidadWorkers; i++) {
            cantidadFilas[i+1] = cantidadInicial;
        }

        for(i=0; i<cantidadExtra; i++) {
            cantidadFilas[i+1] += 1;
        }

    }

    // Difundimos el tamaño del nodo
    MPI_Bcast(&tamanio, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if(globalRank != 0) {
        autoVector = crearVector(tamanio);
    }

    // Le decimos a cada nodo cuantas filas va a recibir
    MPI_Scatter(cantidadFilas, 1, MPI_INT, &cantidadFilasNodo, 1, MPI_INT, 0, MPI_COMM_WORLD);

    // Enviamos las filas a cada nodo
    if (globalRank == 0) {

        int contadorFilas = 0;

        for(i=0; i<cantidadWorkers; i++) {

            for (j=0; j<cantidadFilas[i+1]; j++) {
                MPI_Send(matrizModificada->elementos[contadorFilas], tamanio, MPI_DOUBLE, i+1, 1000, MPI_COMM_WORLD);
                contadorFilas++;
            }

        }

    } else {

        filasNodo = malloc(sizeof(double *) * cantidadFilasNodo);
        resultadoNodo = malloc(sizeof(double) * cantidadFilasNodo);

        for(i=0; i<cantidadFilasNodo; i++) {
            filasNodo[i] = malloc(sizeof(double) * tamanio);
            MPI_Status status;
            MPI_Recv(filasNodo[i], tamanio, MPI_DOUBLE, 0, 1000, MPI_COMM_WORLD, &status);
        }

    }

    do {

        if (globalRank == 0) {

            // Buffereamos el lambda de la iteracion anterior
            lambdaAnterior = lambdaActual;

        }

        // Difundimos el vector de la iteracion actual
        MPI_Bcast(autoVector->elementos, tamanio, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        if(globalRank == 0) {

            int resultadoOffset = 0;

            // ESPERAR EL RESULTADO
            for(i=0; i< cantidadWorkers; i++) {

                if(cantidadFilas[i+1] != 0) {
                    MPI_Status status;
                    MPI_Recv(resultado->elementos + resultadoOffset, cantidadFilas[i+1], MPI_DOUBLE, i+1, 2000+i+1, MPI_COMM_WORLD, &status);
                    resultadoOffset += cantidadFilas[i+1];
                }

            }

        } else {

            int k;

            // MULTIPLICAR
            for (k=0; k<cantidadFilasNodo; k++) {

                resultadoNodo[k] = 0;

                for (j=0; j<tamanio; j++) {
                    resultadoNodo[k] += filasNodo[k][j] * autoVector->elementos[j];
                }

            }

            if (cantidadFilasNodo != 0) {

                // Enviamos los resultados
                MPI_Send(resultadoNodo, cantidadFilasNodo, MPI_DOUBLE, 0, 2000+globalRank, MPI_COMM_WORLD);

            }

        }

        if(globalRank == 0) {

            // Seleccionamos un elemento como pivot (si no se hizo previamente)
            if (pivote<0) {
                pivote = buscarMayorElemento(resultado);
            }

            // Seleccionamos el nuevo autovalor a partir del pivote seleccionado
            lambdaActual = resultado->elementos[pivote];

            // Calculamos el nuevo autovector
            for (i=0;i<tamanio;i++) {
                autoVector->elementos[i] = resultado->elementos[i] / lambdaActual;
            }

            // Calculamos el error actual
            errorActual = fabs((lambdaActual - lambdaAnterior) / lambdaActual);

            // Incrementamos la cantidad de iteraciones en 1
            contadorIteraciones++;

            /*
             * Iteramos hasta que se cumpla una de la siguientes condiciones:
             * 1) La diferencia porcentual entre el autovalor anterior y el actual
             *    sea menor que el error indicado por parámetro.
             * 2) Alcancemos el límite de iteraciones especificado por parámetro.
             */
            condicionSalida = errorActual <= MIN_ERR
                    || contadorIteraciones == MAX_ITER;

        }

        MPI_Bcast(&condicionSalida, 1, MPI_INT, 0, MPI_COMM_WORLD);

    } while (!condicionSalida);

    if(globalRank == 0) {

        // Liberamos la memoria reservada
        eliminarVector(resultado);

        // Seteamos el autovalor calculado en el objeto VV_PROPIO
        vvPropio->autovalor = lambdaActual;

        medioVV = vvPropio;
        medioVV->autovalor = (1 / medioVV->autovalor) + s;

        // Ya no necesitamos la matriz modificada
        eliminarMatriz(matrizModificada);

        // Guardamos el momento final
        instanteFinal = time(NULL);

        // Imprimimos los autovalores y sus autovectores
        imprimir_VV_PROPIO("Mayor autovalor", mayorVV);
        imprimir_VV_PROPIO("Autovalor medio", medioVV);
        imprimir_VV_PROPIO("Menor autovalor", menorVV);

        // Calculamos e imprimimos el tiempo que tardó en ejecutarse el programa
        int tiempoEjecucion = difftime(instanteFinal,instanteInicial);
        printf("\nTiempo de ejecucion: %d segundos\n",tiempoEjecucion);

        // Guardamos los autovalores/autovectores calculados
        almacenar_VV_PROPIOs(menorVV,mayorVV,medioVV);

        // Limpiamos la memoria
        eliminar_VV_PROPIO(mayorVV);
        eliminar_VV_PROPIO(menorVV);
        eliminar_VV_PROPIO(medioVV);

        // Imprimimos la despedida
        printf("\n\n====================================================================\n");
        printf("Fin del programa\n");
        printf("====================================================================\n\n");

    }

    // Finalizamos la librería mpich2
    MPI_Finalize();
    return (EXIT_SUCCESS);
    
}
